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Abstract 

Recently, barrier certificates have been introduced to prove the safety of continuous or hybrid dynamical systems. A barrier 
certificate needs to exhibit some barrier function, which partitions the state space in two subsets: the safe subset in which the 
state can be proved to remain and the complementary subset containing some unsafe region. This approach does not require 
any reachability analysis, but needs the computation of a valid barrier function, which is difficult when considering general 
nonlinear systems and barriers. This paper presents a new approach for the construction of barrier functions for nonlinear 
dynamical systems. The proposed technique searches for the parameters of a parametric barrier function using interval analysis. 
Complex dynamics can be considered without needing any relaxation of the constraints to be satisfied by the barrier function. 
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1 Introduction 

Formal verification aims at proving that a certain behav¬ 
ior or property is fulfilled by a system. Verifying, e.g., the 
safety property for a system consists in ensuring that it 
will never reach a dangerous or an unwanted configura¬ 
tion. Safety verification is usually translated into a reach¬ 
ability analysis problem [4,9,11,29,30]. Starting from an 
initial region, a system must not reach some unsafe re¬ 
gion. Different methods have been considered to address 
this problem. One may explicitly compute the reach¬ 
able region and determine whether the system reaches 
the unsafe region [16]. An alternative idea is to compute 
an invariant for the system, i.e., a region in which the 
system is guaranteed to stay [9]. This paper considers a 
class of invariants determined by barrier functions. 
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A barrier function [23,24] partitions the state space and 
isolates an unsafe region from the part of the state space 
containing the initial region. In [24] polynomial barriers 
are considered for polynomial systems and semi-definite 
programming is used to find satisfying barrier functions. 
Our aim is to extend the class of considered problems to 
non-polynomial systems and to non-polynomial barriers. 
This paper focuses on continuous-time systems. 

The design of a barrier function is formulated as a quan¬ 
tified constraints satisfaction problem (QCSP) [6,25]. 
Interval analysis is then used to find the parameters 
of a barrier function such that the QCSP is satisfied. 
More specifically, the algorithm presented in [18] for ro¬ 
bust controller design is adapted and supplemented with 
some of the pruning schemes found in [7] to solve the 
QCSP associated to the barrier function design. 

The paper is organized as follows. Section 2 introduces 
some related work. Section 3 defines the notion of barrier 
functions and formulates the design of barrier functions 
as a QCSP. Section 4 presents the framework developed 
to solve the QCSP. Design examples are presented in 
Section 5. Section 6 concludes the work. 

In what follows small italic letters x represent real vari¬ 
ables while real vectors x are in bold. Intervals [a;] and in- 
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terval vectors (boxes) [x] are represented between brack¬ 
ets. We denote by M the set of closed intervals over R, 
the set of real numbers. Data structures or sets S are in 
upper-case calligraphic. The derivative of a function x 
with respect to time t is denoted by x. 


2 Related -work 

To prove the safety of a dynamical system, different ap¬ 
proaches have been proposed [4,9,11,15,19,30,31]. One 
way is to explicitly compute an outer approximation of 
the reachable region from the initial region, z. e., the set of 
possible values for the initial state. If it does not intersect 
the unsafe region, then the system is safe. In [4,14,30] the 
reachable region is computed for linear hybrid systems 
for a finite time horizon using geometric representation 
such as polyhedra. The reachable region for non-linear 
systems is computed in [28] using an abstraction of the 
non-linear systems by a linear system expressed in a new 
system of coordinates. The reachability of a non-linear 
system is formulated as an optimization problem in [10]. 
In [8], the Picard-Lindelof operator is combined with 
Taylor models to find the reachable region for non-linear 
hybrid systems. The main downside of the reachabil¬ 
ity approach is the introduction of over-approximations 
during the computations which may lead to difficulties 
to decide whether the system is safe. 

An alternative way to address the safety problem is by 
exhibiting an invariant region in which the system re¬ 
mains. If the invariant does not intersect the unsafe re¬ 
gion then the safety of the system is proved. One way to 
find such an invariant is by using stability properties of 
the considered dynamical system [13] and to search for 
a Lyapunov function. In [22] a sum of squares decom¬ 
position and a semi-definite programming approach are 
employed to find a Lyapunov function for a system with 
polynomial dynamics. A template approach is consid¬ 
ered in [27] to find Lyapunov functions using a branch 
and relax scheme and linear programming to solve the 
induced constraints. A more general idea about invari¬ 
ants is introduced in [24]. Instead of looking for a func¬ 
tion that fulfills some stability conditions, a function is 
searched that separates the initial region from the unsafe 
region. This idea is extended in [16] to search for invari¬ 
ants in conjunctive normal form for hybrid systems. 


3 Formulation 


This section recalls the safety characterization intro¬ 
duced in [24] for continuous-time systems using barrier 
functions. 


3.1 Safety for continuous-time systems 

Consider the autonomous continuous-time perturbed 
dynamical system 

x = /(x,d), (1) 

where x e A C R" is the state vector and d G 2? is a 
constant and bounded disturbance. The set of possible 
initial states at t = 0 is denoted Xq C X. There is some 
unsafe subset C A that shall not be reached by the 
system, whatever XqG Aq at time t = 0 and whatever 
d G 27. We assume that classical hypotheses (see, e.g., 
[5]) on / are satisfied so that (1) has a unique solution 
x(t, xq , d) G A for any given initial value xqG Aq at time 
t = 0 and any d G 27. 

Definition 1 The dynamical system (1) zs safe z/Vxq G 
Xq, Vd G 27 and Vt ^ 0, x(t, xq, d) ^ A„. 

3.2 Barrier certificates 

A way to prove that (1) is safe is by the barrier certificate 
approach introduced in [24]. A barrier is a differentiable 
function i? : A —>■ R that partitions the state space A 
into A_ where 23 (x) ^ 0 and A+ where 23 (x) > 0 such 
that Xq C A_ and A„ C A_|_. Moreover, 23 has to be 
such that Vxo G Aq, Vd G 27, Vi > 0, 23(x(i, Xq, d)) ^ 0. 

Proving that 23(x(i,X q, d)) ^ 0 requires an evaluation 
of the solution of (1) for all xq G Xq and d G 27. Alterna¬ 
tively, [24] provides some sufficient conditions a barrier 
function has to satisfy to prove the safety of a dynamical 
system, see Theorem 1. 

Theorem 1 Consider the dynamical system (1) and the 
sets X, 27, Xq and A„. If there exists a function 23 : A —>■ 
R such that 

VxgAo, R(x)^0, (2) 

Vx gA„, B{x) > 0, (3) 

Vx G A, Vd G 27, 

B(x) = 0 ^ ^ ^g^^\ /(x,d)) < 0, (4) 

then (1) zs safe. 

In (4) (.,.) stands for the dot product in R". In Theo¬ 
rem 1, (2) and (3) ensure that Xq C A_, and A^ C A+, 
while (4) states that if x is on the border between A_ 
and A+ (z.e., 23(x) = 0), then the dynamics / pushes 
the state back in A_ whatever the value of the distur¬ 
bance d. 
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3.3 Parametric barrier functions 

The search for a barrier B is challenging since it is over 
a functional space. As in [24], this paper considers bar¬ 
riers belonging to a family of parametric functions (or 
templates) B (x, p) depending on a parameter vector 
p € 7^ C R"*. Then one may search for some parameter 
value p such that B{x, p) satisfies (2)-(4). 

If there is no p £ 7^ such that 7?(x, p) satisfies (2)-(4), 
this does not mean that the system is not safe: other 
structures of functions i7(x, p) could provide a barrier 
certificate. 

4 Characterization using interval analysis 

This section presents an approach to find a barrier func¬ 
tion that fulfills the constraints of Theorem 1. These 
constraints are first reformulated to cast the design of 
a barrier function as a quantified constraint satisfaction 
problem (QCSP) [25]. 

4-.1 Constraint satisfaction problem 

Assume that there exist some functions go : A” —K and 
g„ : A —K, such that 

Ao = {xe A|go(x) ^0} (5) 

and 

A„ = {xe A|g„(x) <0}. (6) 

Theorem 1 may be reformulated as follows. 

Proposition 2 If 3p gV such that Vx £ A, Vd £ V 

?(x,p,d) =(go(x) >0VB(x,p) ^0) (7) 

A (g„(x) > 0 VB(x, p) > 0) (8) 

A ^S(x,p) 7^ 0 V ^|j^(x,p),/(x,d)^ < 0^ (9) 

holds true, then the dynamical system (1) safe. 

PROOF. The first component of f (x, p, d), 

^0 (x,p) = (go(x) > 0 V R(x,p) ^ 0) (10) 

may be rewritten as 

^0 (x, p) = (go(x) 0 R(x, p) < 0), 

see, e.g., [12]. If ^o(x, p) holds true for some p £ 7^ 
and X £ A, then one has either x £ Aq and 7?(x, p) ^ 
0, or X ^ Aq. In both cases, (2) is satisfied. A similar 


derivation can be made for the second component of 
f (x, p, d) to encode (3), 

C«(x,p) = (g«(x) > 0 V R(x,p) > 0). (11) 

Now, one may rewrite the last component of t (x, p, d), 

Cf>(x,P,d) = ^R(x,p) 7 ^ 0 V ^^(x,p),/(x,d)^ < 0^ 

( 12 ) 

as 

6(x,P,d) = 

^R(x,p)=0 ^ {^||^(x,p),/(x,d)^ < 0^ , (13) 

which corresponds to (4). If 3p £ 7^ such that Vx £ A, 
Vd £ 77, ^ (x, p, d) holds true, then the conditions of 
Theorem 1 are satisfied and (1) is safe. □ 

In [24], (9) is relaxed into 

Sb (x,p,d) = ^^|j^(x,p),/(x,d)^ < 0^ , (14) 

with the consequence of possible elimination of barrier 
functions that would satisfy (9) for some p but not (14). 
Our aim in this paper is to design barrier functions with¬ 
out resorting to this relaxation by considering meth¬ 
ods from interval analysis [17] which allow to consider 
strongly nonlinear dynamics and barrier functions. 


4.2 Solving the constraints 

To find a valid barrier function one needs to find some 
p £ 7^ such that i7(x, p) satisfies the conditions of 
Proposition 2. For that purpose, the Computable Suf¬ 
ficient Conditions-Feasible Point Searcher (CSC-FPS) 
algorithm [18] is adapted. 

In what follows, we assume that A, 77, and V are boxes, 
z.e., A = [x], 77 = [d], and T’ = [p]. CSC-FPS may also 
be applied when A, 77, and V consist of a union of non¬ 
overlapping boxes. 

Consider some function g : R" x R™ —>■ R^' and some box 
[z] £ IR^. CSC-FPS is designed to determine whether 

3p £ [p], Vx £ [x], g (x, p) £ [z] (15) 

and to provide some satisfying p. We extend CSC-FPS to 
handle conjunctions and disjunctions of constraints and 
supplement it with efficient pruning techniques involving 
contractors provided by interval analysis [17]. 
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FPS branches over the parameter search box [p]. 
Branching is performed based on the results provided 
by CSC. For a given box [p]p C [p], CSC returns true 
when it manages to prove that (15) is satisfied for some 
p G [p]g. CSC returns false when it is able to show 
that there is no p S [p]g satisfying (15). CSC returns 
unknown in the other cases. 

In Proposition 2, ^ (x, p, d) consists of the conjunction 
of three terms of the form 

T (x, p, d) = (u(x, p) G ^) V (w(x, p, d) G S). (16) 

For ^o(x, p), defined in (7), 

^ = ]0,+oo] and = ]—oo, 0]; (17) 

for ^„(x, p), defined in (8), 

^ = ]0,+oo] and = ]0,+oo[; (18) 

for ^(,(x,p,d), defined in (9), 

= ]—oo, 0[ U ]0,+oo[ and S = ]—oo, 0[. (19) 


To illustrate the main ideas of CSC-FPS combined with 
contractors, one focuses on the generic QCSP 

3p G [p], Vx G [x] , Vd G [d], r (x, p, d) holds true. 

( 20 ) 

Finding a solution for such QCSP involves three steps: 
validation, reduction of the parameter and state spaces, 
and bisection. 

4-2.1 Validation 

In the validation step, one tries to prove that some vector 
p G [p] is such that Vx G [x], Vd G [d], t (x, p, d) holds 
true. By definition of r (x, p, d), one has to prove that 

3p G [p] ,Vx G [x] ,Vd G [d] , 

(u(x,p) G Vl) V (z;(x,p,d) G S). (21) 

For that purpose, one chooses some arbitrary p G 
[p] and evaluates the set of values u([x],p) = 
{u(x,p) I X G [x]} and u([x] ,p, [d]) = {i;(x,p,d) | 
X G [x],d G [d]} taken by u(x, p) and u(x, p, d) for 
all X G [x] and for all d G [d]. Outer-approximations 
of M ([x], p) and u ([x], p, [d]) are easily obtained using 
inclusion functions provided by interval analysis. 

Definition 3 An inclusion function [f] : IR” —>■ 
for a function f : R" —>■ R^ satisfies for all [x] G IR”, 

/([x]) = {/(x)|xG[x]}C[/]([x]). (22) 


The natural inclusion function is the simplest to obtain: 
all occurrences of the real variables are replaced by their 
interval counterpart and all arithmetic operations are 
evaluated using interval arithmetic. More sophisticated 
inclusion functions such as the centered form, or the 
Taylor inclusion function may also be used, see [17]. 

Using inclusion functions, one may evaluate whether 

[u]([x],p) C A or [?;]([x],p, [d]) C B holds true 

for the various A and B defined in (17), (18), and (19). 

Different choices can be considered for p: one can take 
a random point in [p], the middle, or one of the edges of 
[p]. Here, only the middle of [p] is considered. 

4.2.2 Reduction of the parameter and state spaces 

To facilitate the search for p G [p] one may previously 
eliminate parts of [p] which may be proved not to contain 
any p satisfying (21). The elimination process can be 
done by evaluation or by using contractors [7]. 


4.2.2.1 Evaluation From (21) one deduces that a 
box [p] can be eliminated, le., shown not to contain any 
p satisfying (16), if 

Vp G [p],3x G [x] ,3d G [d] , 

u(x, p) C M A u(x, p, d) C H, (23) 

where A = R\Vl and B = R\H. Using again an inclusion 
function, one may verify whether [m](x, [p]) C A and 
M(x, [p] ,d) C B for some x G [x] and d G [d] and thus 
eliminate the box [p]. 

One may easily verify (7) and (8) since A and B are 
unbounded and may contain [u] (x, [p]) and [u] (x, [p], d). 
For (9) the inclusion is impossible to prove except in 
degenerate cases since A = {0}. 


4.2.2.2 Contractors Consider some function 
g : R” — >• R^ and some set Z C R^. 

Definition 4 A contractor Cc ■ IR” —>■ IR” associated 
to the generic constraint 

c : g(x) G Z (24) 

is a function taking a box [x] as input and returning a 
hoxCc{\y4\) satisfying 

Cc([x])C[x] (25) 
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and 


g{[x])r\Z = g{Cc{[x]))r\Z. 


(26) 


«) 


Cc provides a box containing the solutions of ( 7 (x) £ 
Z included in [x]: (25) ensures that the returned box 
is included in [x] and (26) ensures that no solution of 
g{x) G Z in [x] is lost. 

Consider now two functions gi : M" —>■ and g 2 : 

K" — two sets Z\ C R^^ and Z^ C R^^, and 
the associated constraints ci : 5 i([x]) C Zi and C2 : 
52 ([x]) C Z 2 . Assume that two contractors Cc^ and Cc^ 
are available for ci and C 2 . A contractor CciAc 2 for the 
conjunction c\ A C 2 of c\ and C 2 may be obtained as 

CciAo 2 (M) = ([x]) nCc 2 ([x]), (27) 

or by composition of contractors 

CciAc 2 (M) = ^-02 (f-ci ([x])). (28) 

A contractor CciVc 2 for the disjunction c\ V C 2 of ci and 
C 2 may be obtained as follows 

CciVc2(M) = '^{f'ci([x]) UCc2([x])}, (29) 

see [7], with □{•} the interval hull of a set. 

Using a contractor Cc for (24), one is able to characterize 
some [x] C [x] such that Vx G [x], g(x) ^ Z. 

Proposition 5 Consider a box [x], the elementary con¬ 
straint (24), and the contracted boxCc ([x]) C [x]. Then, 

Vx £ [x] \Cc ([x]) , one has g(x) ^ Z, (30) 

where [x]\Cc([x]) denotes the box [x] deprived from 
Cc ([x]), which is not necessarily a box. 


PROOF. Consider x £ [x]\Cc([x]) and assume that 
g{x) £ Z. Since g{x) G Z and x £ [x], one should have 
X £ Cc([x]), according to (26), which contradicts the 
fact that X £ [x] \Cc ([x]). □ 

Proposition 5 can be used to eliminate [p] or a part of 
[p] for which it is not possible to find any p satisfying 
(21). Consider the constraint 

T : (u(x, p) £ A) V (r;(x, p,d) £ B) (31) 

and a contractor Ct for this constraint. It involves ele¬ 
mentary contractors for the components of the disjunc¬ 
tion in (31), combined as in (29). For the boxes [x], [p], 
and [d] , one gets 

([x]', [p]', [d]') =C.([x], [p], [d]), (32) 



[P]’ 

P 



[d] 


[x]', [d]’ 



Fig. 1. Contractions using Ct', the set for which (31) is satis¬ 
fied is in gray; (a) [p] \ [p]^ yf 0 and the search space for sat¬ 
isfying p can be reduced to [p]'; (6) [x]' / [x] or [d]' / [d], 
it is thus not possible to find some p £ [p] such that (31) is 
satisfied for all x £ [x] and all d £ [d]. 

where [x]^, [p]\ and [d]^ are the contracted boxes. Three 
cases may then be considered. 

(1) If [P] \ [P]' 7 ^ 0: then Vp £ [p] \[p]', Vx £ [x], Vd £ 
[d], 

u(x, p) ^ A A ?;(x, p, d) ^ S, (33) 

and there is no p £ [p] \ [p]^ such that (31) holds 
true for all x £ [x] and for all d £ [d]. Consequently, 
the search space for p can be reduced to [p]^ see 
Figure 1 (a). 

(2) If [x] \ [x]^ 0 then, from Proposition 5, one has 

Vp £[p], Vx £[x] \[x]', Vd £[d], 

u(x, p) ^ A A n(x, p, d) ^ S, (34) 

and there is no p £ [p] such that (31) holds true for 
all X £ [x], see Figure 1 (b). 

(3) If [d] \ [d]' ^ 0, then Vp £ [p], Vx £ [x], Vd £ 
[d] \ [d]'. 


u(x, p) ^ A A ?;(x, p, d) ^ S, (35) 

and there is no p £ [p] such that (31) holds true for 
all d £ [d], see Figure 1 (b). 
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One can reduce the size of the sets for the state x and the 
disturbance d on which (31) has to be verified using the 
contraction on the negation of this constraint. Consider 





























Fig. 2. Contractions using C^; the set for which (36) is satis¬ 
fied is in white; (a) [x]” 7 ^ [x] and/or [d]” 7 ^ [d] and one has 
only to find some suitable p € [p] such that (31) is satisfied 
for all X £ [x] ” and all d G [d]”; ( 6 ) [p]” 7 ^ [p], one may 
choose any p € [p] \[p] ” (for example the value of p indicated 
in red) and (31) will hold true for all (x, d) G ([x] x [d]). 

the negation t of r 

T = (m(x, p) G 34 ) a (f (x, p, d) G B) (36) 

and a contractor C— for this constraint. Assume that after 
applying C- for the boxes [x], [p] , and [d] , one gets 

([x]”, [p]”, [d]”) =C^([x], [p], [d]). (37) 

From Proposition 5, one knows that 

V (x, p, d) G ([x] X [p] X [d]) \ ([x] ” X [p] ” X [d] ”) , 
m(x, p) G A V z;(x, p, d) G S. (38) 

Indeed, if [p] ” = [p] , one can focus on the search for some 
p G [p] satisfying (31) by considering only [x]” x [d]”, 
since for all (x, d) G ([x] x [d]) \ ([x]” x [d]”), (31) is 
satisfied for all p G [p], see Figure 2 (a). 

Now, if [p] ” 7 ^ [p] , then any p G [p] \ [p] ” will satisfy 
(31) for all (x, d) G ([x] x [d]), see Figure 2 (&). 

4-2.3 Bisection 

One is unable to decide whether some p G [p] satisfies 
( 21 ) when 

[u]([x],p) n A 7 ^ 0 and [m]([x],p) ^ A (39) 


or when 

H(M,P, [d])n 6 7^0and[t;]([x],p,[d]) (40) 

This situation occurs in two cases. First, when the se¬ 
lected p does not satisfy (21) for all x G [x] and for all d G 
[d]. Second, when inclusion functions introduce some 
pessimism, i.e., they provide an over-approximation of 
the range of functions over intervals. 

To address both cases, one may perform bisections of 
[x] X [d] and try to verify (21) on the resulting sub-boxes 
for the same p. Bisection allows to isolate subsets of 
[x] X [d] on which one may show that (23) holds true. 
Bisections also reduce pessimism, and may thus facili¬ 
tate the verification of (21). The bisection of [x] x [d] 
is performed within CSC as long as the width of the bi¬ 
sected boxes are larger than some Ex > 0. When CSC is 
unable to prove (21) or (23) and when all bisected boxes 
are smaller than e^, CSC returns unknown. 

FPS performs similar bisections on [p] and stops when 
the width of all bisected boxes are smaller than £p > 0 . 

4-2-4 Composition of constraints 

To prove the safety of the dynamical system (1), Propo¬ 
sition 2 shows that one has to find some p G [p] such 
that Vx G [x], Vd G [d], ^(x, p, d) holds true. Since 
$ (x, p, d) consists of the conjunction of three elemen¬ 
tary constraints of the form ( 20 ), validation requires the 
verification of (21) for the same p considering (7), ( 8 ), 
and (9) simultaneously. Invalidation may be performed 
as soon as one is able to prove that one of the constraints 
(7), ( 8 ), or (9) does not hold true using (23). Contrac¬ 
tion may benefit from the conjunction or disjunction of 
these constraints, as introduced in (28) and (29). 

4-3 CSC-FPS algorithms with contractors 

The CSC-FPS algorithm, presented in [18] is supple¬ 
mented with the contractors introduced in Section 4.2 
to improve its efficiency. 

FPS, described in Algorithm 1, searches for some sat¬ 
isfying p G [p], i.e., some p G [p] such that ^ (x,p,d) 
introduced in Proposition 2 holds true for all x G [x] 
and all d G [d] . This may require to bisect [p] into sub¬ 
boxes stored in a queue Q, which initial content is [p]. 
A sub-box [pjg C [p] is extracted from Q. A reduction 
of [p]q is then performed at Line 1 to eliminate values of 
p G [p]g which cannot be satisfying. To facilitate con¬ 
traction, specific X G [x] are chosen; here only the mid¬ 
point m([x]) of [x] is considered. If [pjg is empty, the 
next box in Q has to be explored. Then CSC is called for 
each constraint to, tu, and 4 to verify whether m (jpjg), 
the midpoint of [pjg, is satisfying. When all calls of CSC 
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return true at Line 1, a barrier function with parameter 
m ([p]q) is found. When one of the calls of CSC returns 
false at Line 1, [p]q is proved not to contain any sat¬ 
isfying p. In all other cases, when w([p]q), the width of 
[p]q, is larger than £p, [p]p is bisected and the resulting 
subboxes are stored in Q. When [pjp is too small, even if 
one was not able to decide whether it constains a satisfy¬ 
ing p, it is not further considered to ensure termination 
of FPS in finite time. The price to be paid in such situa¬ 
tions is the impossibility to conclude that the initial box 
[p] does not contain some satisfying p. This is done by 
setting the flag to false at Line 1. Finally, when Q = 0, 
no satisying p has been found. Whether [p] may however 
contain some satisfying p depends on the value of flag. 


Algorithm 1 FPS 

1 : procedure FPS(^o, ^b, [p],[x], [d]) > ^o, Cu, 6 

from (7), (8) and (9) 

2 : queue Q := [p] 

3: flag true 

4: while Q 7 ^ 0 do 

5: [p]o-~ dequeue(Q) > Reduction of the parameter 

space using (28), (32), (33) 

6 : [p]o := (m([x]), [p]^) t> When [p]q = 0, 

there is no satisfying p € [p]g 
7: if [p]q = 0 then 

8 : continue 

9: end if 

10: ro:= CSC(Co,[p]o, [x],[d]) > Call CSC for each 

constraint 

11: r„:= CSC(C.,[pi;,[x],[d]) 

12: CSC(6,[p];,[x],[d]) 

13: if (ro=true)A(ru=true)A(r(,=true) then 

14: return(true,m([p]p) > m([p]p) is satisfying if 

CSCs hold true 
15: end if 

16: if (ro=f alse)V(ru=false)V(r(,=false) then 

17: continue > no solution in [p]o if one 

constraint does not hold true 
18: end if 

19: if w([p]g) < £p then 

20 : flag:=false > no decision could be made for 

[p]o and thus bisect 

2 ^ • 6ls0 

22: ([p]i,[p]2):=bisect([p]') 

23: enqueue([p]j^) in Q 

24: enqueue([p] 2 ) in Q 

25: end if 

26: end while 

27: if flag=false then 

28: return(unknown, 0 ) > precision reached no 

conclusion can be made for [p] 

29: end if 

30: return(f alse,0) > no valid solution in [p] 

31: end procedure 


CSC, described in Algorithm 2, determines either 
whether r (x, m([p]), d) introduced in (16) holds true 
for all x £ [x] and all d £ [d] or whether there is no 
p £ [p] satisfying (16) for all x £ [x] and all d £ [d]. 
For that purpose, due to the pessimism of inclusion 


functions, it may be necessary to bisect [x] x [d] in sub¬ 
boxes stored in a stack S initialized with [x] x [d]. For 
each subbox [x]p x [d]p C [x] x [d], CSC determines 
at Line 7 whether m([p]Q) is satisfying. CSC tries to 
prove that [p]g does not contain any satisfying p. This 
is done at Line 10, for some (x, d) £ [x]p x [d]p, here 
taken as the midpoint of [x]q x [djg, by an evaluation of 
r. This is done at Line 12 using the result of a contrac¬ 
tor, as described in (34) and (35). When one is not able 
to conclude and provided that w([x]q x [d]g) is larger 
than £x, some parts of [x]p x [d]g for which m([p]g) is 
satisfying are removed at Line 19, before performing a 
bisection and storing the resulting subboxes in S. When 
w([x]p X [d]p) is less than it is no more considered. 
As a consequence, one is not able to determine whether 
m([p]g) is satisfying for all (x, d) £ [x]q x [djp. One may 
still prove that one [p]q does not contain any satisfying 
p considering other subboxes of [x] x [d], but not that 
m([p]Q) is satisfying for all (x, d) £ [x] x [d]. This is 
indicated by setting flag to unknown at Line 17. 

Algorithm 2 CSC 

1: procedure CSC(r,[p]p,[x],[d]) > r is of the form (16) 
2: stack S := [x] x [d] 

3: flag:=true 

4: while S 7 ^ 0 do 

5: [x]q X [d]o :=pop(5') 

6: if M(Mo, m([p]o)) ^ 

'4vM([x](,, m([p]Q), [d]g) C B then 
7: continue > Validation using (21) 

8 : end if 

9: if M (m ([x]p) , [pjp) n A = 0 A 

M (m (Mo) . [p]o- m (fdlo)) n B = 0 then 
10: return(f alse) > Reduction of the parameter 

space using (23) 

11 : end if 

12: (Mo. [p]o. Mo) := Cr(Mo. Mo- Mo)^ Reduction 

of the parameter space applying (32) 

13: if [x]o 7 ^ [x](, V [d]o / [d],, then 

14: return(f alse) > Consequence of (34), (35) 

15: end if 

16: if (w([x]p X [d]p) < Ex) then 

17: flag:=unknown 

18: else 

19: ([x]", [d]") := C-([x]o, m([p](,), [dJJ > Reduc¬ 

tion of the state space using (37), (38) and bisection 
20: ([x]^ X [d]^, [x ]2 X [d] 2 ):=bisect([x]" x [d]") 

21: push([x]j^ X [d]j^) in S 

22: push([x ]2 x [d] 2 )in S 

23: end if 

24: end while 

25: return(flag) 

26: end procedure 


5 Examples 

This section presents experiments for the characteriza¬ 
tion of barrier functions. The considered dynamical sys¬ 
tems are described first before providing numerical re¬ 
sults, comparison of different approaches and discus¬ 
sions. 
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5.1 Dynamical system descriptions 


For the following examples, one provides the dynamics 
of the system, the constraints go and Qu for the definition 
of the sets Xq and X ^, the state space [x], and the para¬ 
metric expression of the barrier function. In all cases, 
the parameter space is chosen as [p] = [—10, lO]'” where 
m is the number of parameters. 

Example 1 Consider the system 


Xi + X2 

^XiX2 — 0.5a:|y 


with (?o(x) = {xi + 1.25)^ -I- {x 2 — 1.25)^ — 0.05, gu(x) = 
{xi + 2.5)2 -t {x 2 - 0.8)2 _ 0.05, and [x] = [-10^,0] x 
[—10^,10^]. The parametric barrier function is = 

PlP2{xo+P3) 

(.xo+P3)^+pI 


-f Xi -|-P4- 


Example 2 Consider the system from [2] 


Xi \ l—Xi + XiX 2 

X2J V “2:2 

with goipc) = {xi — 1.125)2 -|- [x 2 — 0.625)2 — 0.0125, 
gufx.) = {xi — 0.875)2 _|_ _ q. 125)2 — 0.0125, and 

[a;] = [—100,100] x [—100,100]. The parametric barrier 
function used is Bfx.,p) = \ii(jiiXi) — ln{x 2 )+P 2 X 2 +P 3 - 


with gofx.) = {xi - 1)2 -|- (x 2 + 1.5)2 - 0.05, 5 „(x) = 
(xi + 0.6) + (a ;2 — 1)2 — 0.05, and [x] = [—10^, 10^] x 
[—10^,10^]. The parametric barrier function used is 

B(x,p) = (^)%(^)'- 1 . 

Example 6 Consider the Lorenz system from [32] 




^ 10(0:2-Xl) \ 

a;2 

= 

xi (28 — X3) — X 2 



\ X1X2 - |x 3 y 


with go(x) = (a;i + 14.5)2 -|- (x 2 -f 14.5)2 -|- ( 0:3 — 12.5)2 — 
0.25, g„(x) = {xi + 16.5)2 -h (a ;2 -b 14.5)2 -b (13 - 2.5)2 - 
0.25, and [x] = [—20,20] x [—20,0] x [—20,20]. The con¬ 
sidered parametric barrier function is i?(x, p) = pia ;2 -b 

P2X1 +poXz -bp 4 - 

Example 7 Consider the system from [20] 




^ — X\ -b 4X2 — 6 X 3 X 4 ^ 

X 2 


— Xl — X2 + Xg 

X 3 


X 1 X 4 — X 3 -b X 4 X 6 

X 4 


X 1 X 3 + X 3 X 6 - X4 

X 5 


— 2 x| — X 5 -b Xe 

[xej 


\ - 3 X 3 X 4 - Xg - X6 / 


Example 3 Consider the system from [21 [ 


Xl 

yX2, 


X2 

xl-\-x 2 


y'l + {xl+x2y ^ 

with go(x) = a:i2 -b 0:22 — 0.5, gu(x) = {xi — 3.5)2 _|_ 
(a ;2 - 0.5)2 _ Q _5 ^ ^ [-103,103] x [-100,100]. 

A quadratic parametric barrier function is chosen 
B(pc,p) = Pix\ +P 2 x\ +P 0 X 1 X 2 -bp4a;i +P 5 X 2 +Pq. 


with gofx.) = {xi — 3.05)2 _|_ _ 3 . 05^2 _|_ _ 3 . 95^2 _|_ 

(x4-3.05)2-b(a;5-3.05)2-b(x6-3.05)2-0.0001,g„(x) = 
(xi - 7.05)2 - 3.05)2 _ 7 05)2 _ 7 05)2 

(X 5 - 7.05)2 _ 7 05)2 _ 0 0001 ^ and [x] = [0,10] x 

[0,10] X [2,10] X [0,10] X [0,10] X [0,10]. The considered 
parametric barrier function is i?(x, p) = pixl + P 2 X 2 + 
Psxl + Pixl + pox\ + poxl + pr- 

5.2 Experimental conditions and results 


Example 4 Consider the disturbed system from [24] 

Xl\ _ / X 2 \ 

X 2 J \-Xi + ^xf- X 2 J 

with go^x.) = (a:i — 1.5)2 -b 2:22 — 0.25, gu(x) = {xi + 
0.8)2-b(a; 2 -bl)^-0.25, ^ [-100,100] x [-10,10], and 

d€ [0.9,1.1]. The parametric barrier function B {x., p) = 
Pixf + P 2 X 2 + P 3 X 1 X 2 -bp 4 a;i +^ 5 X 2 + Po is considered. 

Example 5 Consider the system with a limit cycle 

(— ( ^2 + (1 - a:f - xl)xi + ln(xf + 1) \ 

\i 2 / \—a;i -b (1 — a:i — X 2 )a ;2 -b ln(a ;2 + 1 )/ 


CSC-FPS, presented in Section 4, has been implemented 
using the IBEX library [3,7]. The selection of candi¬ 
date barrier functions is performed choosing polynomi¬ 
als with increasing degree, except for Examples 1,2, and 
4, where parametric functions taken from [1,2] are con¬ 
sidered. 

For each example, the computing time to get a valid bar¬ 
rier function and the number of bisections of the search 
box [p] are provided. Table 1 summarizes the results for 
the versions of CSC-FPS with and without contractors. 
As in [18], we choose sTx = 10“^ and £p = 10“®. Com¬ 
putations were done on an Intel core 1.7 GHz processor 
with 8 GB of RAM. If after 30 minutes of computations 
no valid barrier function has been found, the search is 
stopped. This is denoted by T.O. (time out) in Table 1. 
















Fig. 3. Results for Example 1. 
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Fig. 4. Results for Example 4 with various values of the 
disturbance d. 



Fig. 5. Results for Example 5. 

Moreover, for Examples 1,4, and 5, graphical represen¬ 
tation of the computed barrier functions are provided. 
In Figures 3, 4, and 5, Xq is in green, X^ is in red, the 
bold line is the barrier function and some trajectories 
starting from Xq are also represented. 

The results in Table 1 show the importance of contrac¬ 
tors which are beneficial in all cases. Thanks to contrac¬ 
tors, valid barrier functions were obtained for all exam¬ 
ples, which is not the case employing the original version 
of CSC-FPS proposed in [18]. In theory, both FPS and 
CSC are of exponential complexity in the dimension m 
of the parameter space and n of the state space. In prac¬ 
tice, contractors allow, on the considered examples, to 
break this complexity and to consider high-dimensional 


Table 1 

Results of CSC-FPS without and with contractors 



Without contr. 

With contr. 

Example 

n 

m 

time 

bisect. 

time 

bisect. 

1 

2 

4 

36s 

4520 

16s 

4553 

2 

2 

3 

T.O. 

/ 

Is 

159 

3 

2 

6 

1133s 

20388 

Is 

6 

4 

2 

6 

253s 

14733 

7s 

435 

5 

2 

4 

T.O. 

/ 

98s 

4072 

6 

3 

4 

167s 

1753 

21 s 

47 

7 

6 

7 

697s 

67600 

Is 

261 


problems. 

The search for barrier functions using the relaxed ver¬ 
sion (14) of the constraint (12) as in [24], has also been 
performed using the version of our algorithms with con¬ 
tractors considering the same parametric barrier func¬ 
tions. We were not able to find a valid barrier function in 
less than 30 minutes. This shows the detrimental effect 
of the relaxation (14) on the search technique. 

Some examples were also addressed using RSolver, which 
is a tool to solve some classes of QCSP [25]. Nevertheless, 
this requires some modifications of the examples, since 
RSolver does not address dynamics or constraints in¬ 
volving divisions [26]. Moreover, the type of constraints 
processed by RSolver, for problems such as (20), allows 
only parametric barrier functions that are linear in the 
parameters. Outer-approximating boxes for the initial 
and the unsafe regions ftb, X^ defined by go and g„ are 
given to Rsolver. As a consequence. Examples 1, 2 and 
3 could not be considered. Only Examples 4, 5, and 6 
were tested by Rsolver, which was able to find a satis¬ 
fying barrier function only for Lorenz in less than Is. 
RSolver was unable to find a solution for Examples 4 and 
5. Rsolver is well designed for problems with paramet¬ 
ric barriers linear in the parameters, but has difficulties 
for non-linear problem and non-convex constraints such 
as (12). 

6 Conclusion 

This paper presents a new method to find paramet¬ 
ric barrier functions for nonlinear continuous-time per¬ 
turbed dynamical systems. The proposed technique has 
no restriction regarding the dynamics nor the form of 
the barrier function. The search for barrier functions 
is formulated as an interval quantified constraint stat- 
isfaction problem. A branch-and-prune algorithm pro¬ 
posed in [18] has been supplemented with contractors to 
address this problem. Contractors are instrumental in 
solving problems with large number of parameters. The 
proposed approach can thus hnd barrier functions for a 
large class of possibly perturbed dynamical systems. 
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Alternative techniques based on RSolver may be signifi¬ 
cantly more efficient for some specific classes of problems 
where the parameters appear linearly in the parametric 
barrier functions. A combination of RSolver and our ap¬ 
proach may be useful to improve the global efficiency of 
barrier function caracterization. 

Future work will be dedicated to the search for the class 
of parametric barrier functions to consider. This may be 
done by exploring a library of candidate barrier func¬ 
tions. In our approach rejection of a candidate function 
occurs mainly after a timeout. Even if contractors aim¬ 
ing at eliminating some parts of the parameter space 
were defined, their efficiency is limited. Better contrac¬ 
tors for that purpose may be very helpful. 

An other future research direction is to extend the pro¬ 
posed method to hybrid systems as done in [24], z.e., to 
consider a set of quantified constraints for each location 
of an hybrid automaton and the constraints associated 
to the transitions. 
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